function [results_vars, results_vectors] = optimtax_hsv(vers, knitro_path, varargin)
% compute outcomes under crp tax system

parser = inputParser;
parser.addParameter('drop_start',0,@isnumeric)
parser.addParameter('drop_stop',0.7,@isnumeric)
parser.addParameter('drop_step',0.01,@isnumeric)
parser.addParameter('compute_compensation',false,@islogical)
parser.addParameter('compensation_init_file',[],@ischar)
parser.addParameter('welfare_target',[],@isnumeric)
parser.addParameter('welfare_targets_file',[],@ischar)
parser.addParameter('wage_and_price_file',[],@ischar)
parser.addParameter('suffix','',@ischar)

% init_file is used for fixing tax progressivity in HSV tax schedule

% to make results comparable with those obtained by GPOPS, load
% wage_and_price_file which makes sure that the same wage bounds
% and robot prices are used as in GPOPS

parser.parse(varargin{:})
drop_start = parser.Results.drop_start
drop_stop = parser.Results.drop_stop
drop_step = parser.Results.drop_step
compute_compensation = parser.Results.compute_compensation
welfare_range = parser.Results.welfare_target
welfare_targets_file = parser.Results.welfare_targets_file
compensation_init_file = parser.Results.compensation_init_file
wage_and_price_file = parser.Results.wage_and_price_file
suffix = parser.Results.suffix

if isempty(wage_and_price_file)
    wage_and_price_file = welfare_targets_file;
end

if ~strcmp(suffix,'')
    suffix = ['_',suffix];
end

file_suffix = [vers,suffix];

addpath(knitro_path);

addpath('../')
addpath('../Densities')
addpath('../UtilityFunctions')
addpath('../ProductionFunctions')
addpath('../Economies/')

calibration_version = '20191107_nochoice';
prefix = calibration_version;
prefix_pf = prefix;
prefix_skill = prefix;

dir_in = '../../../data/generated/matlab/Calibration/';

if compute_compensation
    if contains(welfare_targets_file, 'hsv')
        dir_out = ['../../../data/generated/matlab/OptimTax/' ...
            'compensation_hsv/'];
    else
        dir_out = ['../../../data/generated/matlab/OptimTax/' ...
            'compensation_gpops/'];
    end
else
    dir_out = '../../../data/generated/matlab/OptimTax/hsv/';
end

inequ_par = tools.csv2struct('../../../data/generated/matlab/Calibration/inequ_par.csv');
constants = tools.csv2struct('../../../data/raw/constants.csv');

% load calibrated parameters for skill distribution
par_dist = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_sol.csv']);
mom_dist = tools.csv2struct([dir_in,prefix_skill,'_abilityDistThreeFixC_mom.csv']);

DensSkill = abilityDistFixed(par_dist);

DensPart = ParticipationNormal(par_dist);

% utility function object
par_util.epsilon = constants.epsilon;
Util = UtilityQuasiLinear(par_util);

% combine in SkillAndParticipation object
par_dist.dens_skill = DensSkill;
par_dist.dens_part = DensPart;
par_dist.util = Util;

par_dist.transfer = par_dist.transfer;
par_dist.num_nodes_phi = 100;

par_dist.hsv_tax = true

Dist = SkillAndParticipation(par_dist);

% Production function
par_prod_fun =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_sol.csv']);
mom_prod_fun =  tools.csv2struct([dir_in,prefix_pf,'_KrusellRobotsEach_mom.csv']);

par_prod_fun.parse_inputs = true;
par_prod_fun.KrusellRobotsEach_use_properties = true;
ProdFun = KrusellRobotsEach(par_prod_fun);

% Economy object
par_econ.varnames = {'L_M','L_R','L_C', ...
    'K_B_M','K_B_R','K_B_C', 'K_E','K_S',...
    'tau_B','tau_E','tau_S','transfer','tau_hsv','lambda_hsv','compensation'};
par_econ.pf = ProdFun;
par_econ.util = Util;

par_econ.lb = 1e-16 .* ones(1,8);
par_econ.ub = 1e16 .* ones(1,8);

% robot prices and taxes
inp.tau_B = par_prod_fun.tau_B;
inp.tau_E = par_prod_fun.tau_E;
inp.tau_S = par_prod_fun.tau_S;

par_econ.q_B = mom_prod_fun.q_B;
par_econ.q_E = mom_prod_fun.q_E;
par_econ.q_S = mom_prod_fun.q_S;

q_B_0 = par_econ.q_B;

par_econ.dist = Dist;

par_econ.solver = 'knitro_eqsolve';

par_econ.phi_lb = constants.phi_lb;
par_econ.r = inequ_par.inequ_par;
par_econ.par_dist = par_dist;

par_econ.w_lb = constants.w0_gpops;
par_econ.w_ub = constants.wf_gpops;

econ = KrusellWithRobotsEachEconomy(par_econ);
econ.set_num_nodes(100,50,20,50);

% compute mass for normalization
integrand = @(w) econ.dist.dens_skill.f(w,par_dist.Y_M_0,par_dist.Y_R_0,par_dist.Y_C_0,par_dist);
% integrand = @(w)
% auxdata.econ.dist.dens_skill.f(w,6.3720,11.7917,19.6443,par_dist);
w0 = par_econ.w_lb;
wf = par_econ.w_ub;
mass = integral(integrand,w0,wf);

econ.dist.dens_skill.mass = mass;

integrand = @(w) econ.dist.dens_skill.f(w,par_dist.Y_M_0,par_dist.Y_R_0,par_dist.Y_C_0,par_dist);
mass_normalized = integral(integrand,w0,wf)

inp.L_M = par_prod_fun.L_M_0 + 1e-3;
inp.L_R = par_prod_fun.L_R_0 + 1e-3;
inp.L_C = par_prod_fun.L_C_0 + 1e-3;
inp.K_B_M = par_prod_fun.K_B_M_0 + 1e-3;
inp.K_B_R = par_prod_fun.K_B_R_0 + 1e-3;
inp.K_B_C = par_prod_fun.K_B_C_0 + 1e-3;
inp.K_E = par_prod_fun.K_E_0;
inp.K_S = par_prod_fun.K_S_0;
inp.transfer = par_dist.transfer;
inp.tau_hsv = par_dist.tau_hsv;
inp.lambda_hsv = par_dist.lambda_hsv;
inp.compensation = 0;

econ.revenue_requirement = mom_prod_fun.revenue_requirement;

switch vers
    case 'adjust_tau_hsv'
        lb.L_M = 1e-6;
        lb.L_R = 1e-6;
        lb.L_C = 1e-6;
        lb.K_B_M = 1e-12;
        lb.K_B_R = 1e-12;
        lb.K_B_C = 1e-12;
        lb.K_E = 1e-6;
        lb.K_S = 1e-6;
        lb.tau_B = inp.tau_B;
        lb.tau_E = inp.tau_E;
        lb.tau_S = inp.tau_S;
        lb.transfer = 0;
        lb.tau_hsv = 0.05;
        lb.lambda_hsv = 1;
        
        ub.L_M = 1e6;
        ub.L_R = 1e6;
        ub.L_C = 1e6;
        ub.K_B_M = 1e8;
        ub.K_B_R = 1e8;
        ub.K_B_C = 1e8;
        ub.K_E = 1e6;
        ub.K_S = 1e6;
        ub.tau_B = inp.tau_B;
        ub.tau_E = inp.tau_E;
        ub.tau_S = inp.tau_S;
        ub.transfer = 1e6;
        ub.tau_hsv = 1;
        ub.lambda_hsv = 1e5;
        
    case 'baseline'
        % only adjust lambda_hsv and transfer_hsv
        tmp = readtable(compensation_init_file);
        init = table2struct(tmp,'ToScalar',true);
        lb.tau_hsv = init.tau_hsv(1);
        ub.tau_hsv = init.tau_hsv(1);
        
        lb.L_M = 1e-6;
        lb.L_R = 1e-6;
        lb.L_C = 1e-6;
        lb.K_B_M = 1e-12;
        lb.K_B_R = 1e-12;
        lb.K_B_C = 1e-12;
        lb.K_E = 1e-6;
        lb.K_S = 1e-6;
        lb.tau_B = inp.tau_B;
        lb.tau_E = inp.tau_E;
        lb.tau_S = inp.tau_S;
        lb.transfer = 0;
        lb.lambda_hsv = 1;
        
        ub.L_M = 1e6;
        ub.L_R = 1e6;
        ub.L_C = 1e6;
        ub.K_B_M = 1e8;
        ub.K_B_R = 1e8;
        ub.K_B_C = 1e8;
        ub.K_E = 1e6;
        ub.K_S = 1e6;
        ub.tau_B = inp.tau_B;
        ub.tau_E = inp.tau_E;
        ub.tau_S = inp.tau_S;
        ub.transfer = 1e6;
        ub.lambda_hsv = 1e5;
  
  case 'adjust_robot_tax'
        tmp = readtable(compensation_init_file);
        init = table2struct(tmp,'ToScalar',true);
        lb.tau_hsv = init.tau_hsv(1);
        ub.tau_hsv = init.tau_hsv(1);
        
        lb.L_M = 1e-6;
        lb.L_R = 1e-6;
        lb.L_C = 1e-6;
        lb.K_B_M = 1e-12;
        lb.K_B_R = 1e-12;
        lb.K_B_C = 1e-12;
        lb.K_E = 1e-6;
        lb.K_S = 1e-6;
        lb.tau_B = -1;
        lb.tau_E = inp.tau_E;
        lb.tau_S = inp.tau_S;
        lb.transfer = 0;
        lb.lambda_hsv = 1;
        
        ub.L_M = 1e6;
        ub.L_R = 1e6;
        ub.L_C = 1e6;
        ub.K_B_M = 1e8;
        ub.K_B_R = 1e8;
        ub.K_B_C = 1e8;
        ub.K_E = 1e6;
        ub.K_S = 1e6;
        ub.tau_B = 10;
        ub.tau_E = inp.tau_E;
        ub.tau_S = inp.tau_S;
        ub.transfer = 1e6;
        ub.lambda_hsv = 1e5;        

  case 'adjust_taus'
        tmp = readtable(compensation_init_file);
        init = table2struct(tmp,'ToScalar',true);
        lb.tau_hsv = init.tau_hsv(1);
        ub.tau_hsv = init.tau_hsv(1);
        
        lb.L_M = 1e-6;
        lb.L_R = 1e-6;
        lb.L_C = 1e-6;
        lb.K_B_M = 1e-12;
        lb.K_B_R = 1e-12;
        lb.K_B_C = 1e-12;
        lb.K_E = 1e-6;
        lb.K_S = 1e-6;
        lb.tau_B = -1;
        lb.tau_E = -1;
        lb.tau_S = -1;
        lb.transfer = 0;
        lb.lambda_hsv = 1;
        
        ub.L_M = 1e6;
        ub.L_R = 1e6;
        ub.L_C = 1e6;
        ub.K_B_M = 1e8;
        ub.K_B_R = 1e8;
        ub.K_B_C = 1e8;
        ub.K_E = 1e6;
        ub.K_S = 1e6;
        ub.tau_B = 10;
        ub.tau_E = 10;
        ub.tau_S = 10;
        ub.transfer = 1e6;
        ub.lambda_hsv = 1e5;             
        
    otherwise
        error('non-existent version')
end

if compute_compensation
    lb.compensation = -1e6;
    ub.compensation = 100;
else
    lb.compensation = 0;
    ub.compensation = 0;
end

if compute_compensation && ~isempty(welfare_targets_file)
    tmp = readtable(welfare_targets_file);
    targets = table2struct(tmp,'ToScalar',true);
    welfare_range = targets.welfare(:)';
    if strcmp(vers,'baseline')
        lb.compensation = -1e4;
    end
end

% read w's and robot price range from file
if isempty(wage_and_price_file)
    w0_range = par_econ.w_lb;
    wf_range = par_econ.w_ub;
    price_range = econ.q_B;
    drop_range = 0;
else
    wage_and_price_file
    tmp = readtable(wage_and_price_file);
    range_file = table2struct(tmp,'ToScalar',true);
    price_range = range_file.q_B(:)';
    w0_range = range_file.w0;
    wf_range = range_file.wf;
    try
        drop_range = range_file.drop(:)';
    catch
        drop_range = range_file.pricedrop(:)';
    end
end

if compute_compensation
    prefix = [prefix,'_compensation'];
end

dir_out
prefix
file_suffix

file_vars = [dir_out,prefix,'_',file_suffix,'.csv']
file_percentiles = [dir_out,prefix,'_',file_suffix,'_percentiles.csv'];

% clear results files
if exist(file_vars, 'file')==2
    delete(file_vars);
end

if exist(file_percentiles, 'file')==2
    delete(file_percentiles);
end

results = simulate(econ, inp, lb, ub, q_B_0, price_range,welfare_range,prefix, ['simulation_optimtax_',vers],dir_out,par_dist);

    function [results, results_percentiles] = simulate(econ,inp,lb,ub, q_B_0,price_range,welfare_range,prefix, suffix, dir_out, par)
        results_vars = table();
        results_vectors = table();
        
        sol = inp;
        stats = struct();
        percentiles = struct();
        results = table();
        results_percentiles = table();
        
        % compute number of iterations based on drop_stop
        no_iterations = sum(drop_range <= drop_stop);
        
        for i = 1:no_iterations
            % change price of robots and recompute
            price = price_range(i);
            econ.q_B = price;
            
            drop = 1 - price/q_B_0;
            
            fprintf('\n=================================================\n')
            fprintf('%s, price is %f\n', suffix, econ.q_B)
            fprintf('=================================================\n')
            
            tic
            if compute_compensation
                
                welfare_target = welfare_range(i);
                fprintf('welfare target =  %f\n', welfare_target)
                econ.w_lb = w0_range(i);
                econ.w_ub = wf_range(i);
                [sol,fval,exitflag,output] = econ.compute_compensation(sol,lb, ...
                    ub,welfare_target);
            else
                econ.w_lb = w0_range(i);
                econ.w_ub = wf_range(i);
                [sol,fval,exitflag,output] = econ.compute_optimtax(sol,lb, ...
                    ub);
            end
            toc
            
            [stats, percentiles] = compute_stats(econ, sol, par);
            stats.pricedrop = drop;
            stats.exitflag = exitflag;
            stats.fval = fval;
            
            try
                % load results file, provided it exists
                results_vars = struct2table(tools.csv2struct(file_vars));
                results_vectors = ...
                    struct2table(tools.csv2struct(file_percentiles));
            end
            
            results_vars = [results_vars; struct2table(stats)];
            results_vectors = [results_vectors; struct2table(percentiles)]; ...
                
        results_vars_struct = table2struct(results_vars);
        results_vectors_struct = table2struct(results_vectors);
        
        tools.struct2csv(results_vars_struct,file_vars)
        tools.struct2csv(results_vectors_struct,file_percentiles)
        end
        
    end

    function [stats, percentiles_strct] = compute_stats(econ, sol, par)
        stats = sol;
        stats.w0 = econ.w_lb;
        stats.wf = econ.w_ub;
        stats.K_B = sol.K_B_M + sol.K_B_R + sol.K_B_C;
        stats.q_B = econ.q_B;
        stats.q_E = econ.q_E;
        stats.q_S = econ.q_S;
        
        stats.tau_return_B = 1 - 1/(1+stats.tau_B);
        stats.tau_return_E = 1 - 1/(1+stats.tau_E);
        stats.tau_return_S = 1 - 1/(1+stats.tau_S);
        
        [stats.Y_M,stats.Y_R,stats.Y_C,stats.Y_B_M,stats.Y_B_R,stats.Y_B_C,stats.Y_E,stats.Y_S]=econ.factor_prices(sol,par);
        
        w_lb = 1e-6;
        w_ub = 1e4;
        
        stats.part = econ.dist.mass_participants(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,par);
        stats.emp_M = econ.dist.emp_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        stats.emp_R = econ.dist.emp_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        stats.emp_C = econ.dist.emp_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        
        stats.avg_wage_M = econ.dist.average_wage_M(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        stats.avg_wage_R = econ.dist.average_wage_R(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        stats.avg_wage_C = econ.dist.average_wage_C(stats.Y_M,stats.Y_R,stats.Y_C,w_lb,w_ub,stats.part,par);
        
        % welfare
        derived_vars = econ.compute_vars_from_input(sol);
        stats.welfare_part_M = econ.welfare_participants_M(derived_vars);
        stats.welfare_part_R = econ.welfare_participants_R(derived_vars);
        stats.welfare_part_C = econ.welfare_participants_C(derived_vars);
        
        stats.welfare_non_part = econ.welfare_non_participants(derived_vars);
        stats.welfare_part = econ.welfare_participants(derived_vars);
        stats.welfare = econ.welfare(derived_vars);
        
        % tax revenues
        stats.tax_rev_inc = derived_vars.tax_rev_inc;
        stats.exp_transfer = derived_vars.exp_transfer;
        stats.tax_rev_B = derived_vars.tax_rev_B;
        stats.tax_rev_E = derived_vars.tax_rev_E;
        stats.tax_rev_S = derived_vars.tax_rev_S;
        stats.tax_rev_total = derived_vars.tax_rev_total;
        stats.excess_revenue = derived_vars.excess_revenue;
        stats.budget_constraint = derived_vars.budget_constraint;
        
        stats.Y = derived_vars.Y;
        
        % percentiles
        handle_h = @(w) 1/stats.part .* econ.dist.h(w,stats.Y_M, ...
            stats.Y_R,stats.Y_C,par);
        pct = 0.1:0.1:0.9;
        [percentiles, nodes, weights] = ...
            tools.Integration.compute_percentiles(handle_h,pct,1,200,100);
        
        stats.part_M = econ.part_M(sol,par);
        stats.part_R = econ.part_R(sol,par);
        stats.part_C = econ.part_C(sol,par);
        
        % population shares
        integrand = @(w) econ.dist.dens_skill.f_M(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
        stats.pop_M = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);
        
        integrand = @(w) econ.dist.dens_skill.f_R(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
        stats.pop_R = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);
        
        integrand = @(w) econ.dist.dens_skill.f_C(w,stats.Y_M,stats.Y_R,stats.Y_C,par);
        stats.pop_C = tools.Integration.integral_legendre(integrand,w_lb,w_ub,econ.num_nodes_f,econ.nodes_f,econ.weights_f);
        
        percentiles_strct.percentile = pct';
        percentiles_strct.value = percentiles;
        percentiles_strct.q_B = econ.q_B * ones(size(pct'));
        percentiles_strct.q_E = econ.q_E * ones(size(pct'));
        percentiles_strct.q_S = econ.q_S * ones(size(pct'));
    end
end